Thalamocortical Connections Mature along a Sensorimotor-Association Axis of Cortical Developmental Heterochronicty

Read in Data for Developmental Characterization

Glasser parcel names and mesulam assignments

#Glasser regions with corresponding labels, tract names, and mesulam assignments
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -

#Glasser regions assigned to mesulam zones
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv")
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)
glasser.assignments <- glasser.assignments %>% select(label, tract, orig_parcelname, mesulam_assignment)

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per bin

S.A.axis %>% group_by(SA.axis.bin) %>% do(num.regions = length(.$SA.axis), mean.rank = mean(.$SA.axis)) %>% unnest(cols=c("num.regions", "mean.rank"))
## # A tibble: 5 × 3
##   SA.axis.bin num.regions mean.rank
##         <dbl>       <int>     <dbl>
## 1           1          72      36.5
## 2           2          72     108. 
## 3           3          72     180. 
## 4           4          72     252. 
## 5           5          72     324.

NeuroSynth cognitive atlas term maps

#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) #list of terms (cue willy wonka line "dont you want to know our names?" "cant imagine why it wouldnt matter")
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rh

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Dataset-specific diffusion measure spreadsheets

FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #generated by sample_construction/PNC/tractmeasures_dfs_PNC.R

FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv") #generated by /sample_construction/HCPD/tractmeasures_dfs_HCPD.R

Dataset-specific tract lists

tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

Developmental statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir

#list all FA development files in development results directory for accessing. Files were generated by gam_functions/fit_ageGams.R
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F) 

filenames <- c()
#generate variable names to assign files data to 
for(name in files){
  Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
  filenames <- append(filenames, Rname) #save filenames into a character vector 
}

#read in files and assign to variables
for(i in 1:length(filenames)){
  
  Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
  
  if(grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, glasser.assignments, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, S.A.axis, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(grepl("sexeffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, glasser.assignments, by = "tract")
    assign(Rfilename, x) 
  }
  
}
rm(x)

A Spectrum of Thalamocortical Developmental Trajectories

Cortex-wide thalamic connectivity developmental profiles

Number of significant developmental effects

#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

PNC

sig.effects(measure = "FA", atlas = "glasser", dataset = "pnc")
## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes in FA"

HCPD

sig.effects(measure = "FA", atlas = "glasser", dataset = "hcpd")
## [1] "In hcpd, 180 thalamocortical connections (78.26 percent) show significant age-related changes in FA"

Number of significant age by sex interactions

#Function to calculate the number and percent of thalamocortical connections showing a significant age*sex interaction
sig.effects.sex <- function(measure, atlas, dataset){
  sexeffects.df <- get(sprintf("sexeffects_%s_%s_%s", measure, atlas, dataset)) 
  sexeffects.df$significant <- p.adjust(sexeffects.df$GAM.int.pvalue, method = c("fdr")) < 0.05 #fdr-corrected 
  sigeffect.totaln <- sum(sexeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(sexeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age by sex interactions", dataset, sigeffect.totaln, sigeffect.percent))
}

PNC

sig.effects.sex(measure = "FA", atlas = "glasser", dataset = "pnc")
## [1] "In pnc, 0 thalamocortical connections (0 percent) show significant age by sex interactions"

HCPD

sig.effects.sex(measure = "FA", atlas = "glasser", dataset = "hcpd")
## [1] "In hcpd, 4 thalamocortical connections (1.74 percent) show significant age by sex interactions"

Thalamocortical connection GAM smooth functions

PNC

smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)

HCPD

smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")

ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  scale_y_continuous(limits = c(-0.02, 0.0138)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)

Region-specific GAM smooths

#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
  
  #Read in participant-level data and model-level fitted values and derivatives 
  participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
  fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA

  #Connection spline plot with participant-level data
  trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
    geom_point(color = alpha(display_color, 0.5), size = .1) +
    geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
    geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted,  ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
    labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
    theme_classic() +
    scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    scale_y_continuous(breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  #Significant derivatives (rate of age-related change) plot 
  derivatives.plot <- ggplot(data = derivatives.df) + 
    geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
    scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color,  na.value = "white") +
    labs(x="\nAge", fill = "Derivative") + 
    theme_classic() +
    scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    theme(
        legend.position = "none",
        axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size = 6, family = "Arial", color = c("black"))) 
  
  allplots <- list(trajectory.plot, derivatives.plot)
  tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)

  return(tract.plot) 
}

Primary motor cortex (4)

Left 4: S-A axis rank = 16, mesulam assignment = idiotypic

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               16.06784
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               15.45184

Lateral parietal (PF)

Left PF: S-A axis rank = 248, mesulam assignment = heteromodal

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1                18.4531
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               16.49456

Lateral prefrontal (8BL)

Left 8BL: S-A axis rank = 342, mesulam assignment = heteromodal

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               21.21106
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               21.91667

Maturational patterns are highly reproducible

Correlation of thalamocortical connection development effect sizes between datasets

Pearson’s R

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hcpd
## t = 16.04, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6677655 0.7902850
## sample estimates:
##       cor 
## 0.7349671

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0

Correlation plot

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.06, 0.15)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)

Correlation of thalamocortical connection age of maturation between datasets

Pearson’s R

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hcpd
## t = 8.6088, df = 153, p-value = 8.379e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4545288 0.6686745
## sample estimates:
##       cor 
## 0.5712442

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 2e-04

Correlation plot

ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
  geom_point(color = c("#9c3a80"), size = 0.4) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
  scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Agematuration_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.94, height = 1.58)

Correlation of developmental magnitude and age of maturation within dataset

PNC

Spearman’s rho

cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$GAM.smooth.partialR2 and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 234341, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8242399

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0

Correlation plot

ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)

HCPD

cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$GAM.smooth.partialR2 and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 137678, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8161002

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0

Early and Late Maturing Thalamocortical Connections

Early and late maturing areas

#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late

matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey

ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  side  region label smooth.increase.offset maturation_bracket
##   <chr> <chr> <chr> <chr> <chr>  <chr>                  <dbl> <chr>             
## 1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L…                   16.9 mid

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Maturationalbrackets_PNC.pdf", device = "pdf", dpi = 500, width = 4.4, height = 4.4)

NeuroSynth decoding of age of maturation maps

Compute cognitive term-development map spatial correlations

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

Developmental timing decoding by Neurosynth cognitive terms

PNC

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))
devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.pnc
##                   term term.correlation
## 1           expectancy        0.4099804
## 2           monitoring        0.3968318
## 3    cognitive_control        0.3822346
## 4            reasoning        0.3317204
## 5  response_inhibition        0.3131374
## 6             strategy        0.3119431
## 7     memory_retrieval        0.2883541
## 8              thought        0.2863203
## 9               memory        0.2799527
## 10              recall        0.2692535
## 11  object_recognition       -0.2919385
## 12   visual_perception       -0.2896572
## 13          morphology       -0.2747329
## 14       consolidation       -0.2571732
## 15                gaze       -0.2507995
## 16           expertise       -0.2367694
## 17        coordination       -0.2367162
## 18          perception       -0.2249922
## 19   facial_expression       -0.2183252
## 20           induction       -0.2049576

HCPD

gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")

devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
  arrange(desc(term.correlation))
devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))

devpattern.neurosynth.hcpd
##                      term term.correlation
## 1              expectancy        0.3933764
## 2  sentence_comprehension        0.3887129
## 3                 thought        0.3767177
## 4       cognitive_control        0.3477497
## 5     response_inhibition        0.3345571
## 6              monitoring        0.3258089
## 7                  recall        0.3033990
## 8      emotion_regulation        0.2986220
## 9                language        0.2927547
## 10              reasoning        0.2664190
## 11           multisensory       -0.3268444
## 12      visual_perception       -0.3132163
## 13     object_recognition       -0.3088589
## 14              induction       -0.2857823
## 15           localization       -0.2604849
## 16         mental_imagery       -0.2576200
## 17             adaptation       -0.2537276
## 18                   gaze       -0.2474086
## 19               fixation       -0.2473204
## 20             navigation       -0.2404549

Overlapping terms across datasets

merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)
##                   term
## 1           expectancy
## 2           monitoring
## 3    cognitive_control
## 4            reasoning
## 5  response_inhibition
## 6              thought
## 7               recall
## 8            induction
## 9                 gaze
## 10   visual_perception
## 11  object_recognition
ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 2.2, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
  scale_y_discrete(
  labels = c("expectancy"=expression(bold(expectancy)),
  "monitoring"=expression(bold(monitoring)),
  "cognitive_control"=expression(bold(cognitive_control)),
  "reasoning"=expression(bold(reasoning)),
  "response_inhibition"=expression(bold(response_inhibition)),
  "thought"=expression(bold(thought)),
  "recall"=expression(bold(recall)),
  "induction"=expression(bold(induction)),
  "object_recognition"=expression(bold(object_recognition)),
  "visual_perception"=expression(bold(visual_perception)),
  "gaze"=expression(bold(gaze)),
  "induction"=expression(bold(induction), parse=TRUE))) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_line(color = c("gray90")),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)

Probability of obtaining 11/20 overlapping words given the set of 123 terms

#Function to compute the number of overlapping terms when drawing term sets of size 20 at random for PNC and HCPD nulls
compute_term_overlap <- function(){
term.subset.pnc <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
term.subset.hcpd <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
num.overlap <- length(intersect(term.subset.pnc, term.subset.hcpd))
return(num.overlap)}

#Null distribution of overlapping term counts based on 10,000 random draws
set.seed(1)
overlap.number <- replicate(10000, compute_term_overlap()) %>% as.data.frame %>% set_names("overlaps")

sprintf("In %s out of 10,000 runs, the term overlap number was equal to or greater than 11", sum(overlap.number$overlaps >= 11))
## [1] "In 0 out of 10,000 runs, the term overlap number was equal to or greater than 11"
sprintf("The permutation-based p value for these overlap counts is %s", sum(overlap.number$overlaps >= 11) / 10000)
## [1] "The permutation-based p value for these overlap counts is 0"
summary(overlap.number$overlaps)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   3.000   3.257   4.000  10.000
ggplot(overlap.number, aes(x = overlaps)) + 
  geom_histogram(bins = 12, fill = c( "#9c3a80"), color = c("white")) +
  theme_classic() +
  scale_x_continuous(limits = c(-1, 11), breaks = c(0, 2, 4, 6, 8, 10)) +
  labs(x = "\nNeurosynth term overlap", y = "Count\n") +
  geom_vline(xintercept = 11, linewidth = .1) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .1), 
  axis.ticks = element_line(linewidth = .1))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_overlap_null.pdf", device = "pdf", dpi = 500, width = 1.3, height = 1.1)

Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis

gameffects_FA_glasser_pnc <- merge(gameffects_FA_glasser_pnc, (S.A.axis %>% select(tract, SA.axis)), by = "tract")
gameffects_FA_glasser_hcpd <- merge(gameffects_FA_glasser_hcpd, (S.A.axis %>% select(tract, SA.axis)), by = "tract")

Age-resolved developmental change plots

PNC

derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)

HCPD

derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

Thalamocortical connections mature at progressively older ages across the S-A axis

PNC

Spearman’s rho

cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$SA.axis and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 683025, p-value = 2.394e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4877183

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.00075

Correlation plot

ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

Brain map

Thalamocortical Age of Maturation

agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey

ggplot() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  17.5, na.value = "gray93", limits = c(16, 19), oob = squish) + 
  new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) 

Sensorimotor-Association Axis

S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)

ggplot() + 
  geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  180, na.value = "gray93") + 
   new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) 

HCPD

Spearman’s rho

cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$SA.axis and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 367590, p-value = 2.934e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5090025

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.00155

Correlation plot

ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_hcpd.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

Thalamocortical maturation along anatomical axes

#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180]  <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)
#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
  df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
  df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes,  by = c("label", "tract", "orig_parcelname"))
  
  #S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - development correlation
  anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

PNC

Anterior-posterior axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.3101209          0.0621 0.4877183 0.0001552829

Medial-lateral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")
##    axis   axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni 0.03024772         0.43765 0.4877183 3.832011e-08

Dorsal-ventral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.1173292         0.26005 0.4877183 0.0002034701
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4877183)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3101209)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1173292)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.03024772)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "D-V", "M-L"))

ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) + 
  geom_bar(stat='identity') + 
  labs(x="") +
  labs(y="\nAge of Maturation Correlation") +
  theme_classic()+
  scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_pnc.pdf", device = "pdf", dpi = 500, width = 1.35, height = 1)

HCPD

Anterior-posterior axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")
##    axis axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.441284         0.02195 0.5090025     0.124339

Medial-lateral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni 0.1124713          0.3695 0.5090025 5.601018e-07

Dorsal-ventral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.0814776         0.38565 0.5090025 1.390815e-05

Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity

Plasticity marker development data

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)
#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)
#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full

library(ciftiTools)

#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels

#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version 
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released 
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)

#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")

#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)

detach("package:ciftiTools", unload=TRUE)
EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)
df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

Fluctuation amplitude decrease onset map

plasticity.dev$agedecline.protracted <- plasticity.dev$amplitude.agedecline
plasticity.dev$agedecline.protracted[plasticity.dev$agedecline.protracted == "NaN"] <- 23 #lets plot regions with most protracted functional expression of plasticity (amplitude never declines!)
snr.exclude <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/SNRmask_glasser360.csv") %>% filter(SNR.mask == 0) #but not those we excluded for low SNR
plasticity.dev$agedecline.protracted[plasticity.dev$label %in% snr.exclude$label] <- NA

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = agedecline.protracted, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(10, 17), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.7, height = 5.4) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

T1/T2 annualized rate of increase map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 1), alpha("#323280", 0.8), alpha("#323280", 0.6), alpha("#323280", 0.4), alpha("#323280", 0.2)), limits = c(0.012, 0.016), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.4, height = 5.13) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

E:I ratio magnitude of developmental decline map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(-0.04, -0.035), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5, height = 5.08) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts

PNC

plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$EI.ageslope
## S = 733425, p-value = 2.329e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4499176

Spatial permutation (spin) based p-value

EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$myelin.annualroc
## S = 1929655, p-value = 3.141e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4472777

Spatial permutation (spin) based p-value

myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$amplitude.agedecline
## S = 712074, p-value = 3.088e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3028346

Spatial permutation (spin) based p-value

amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.pnc, myelin.thalamus.pnc, amplitude.thalamus.pnc)
plasticity.ps
## [1] 0.00020 0.00090 0.02965
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.00060 0.00135 0.02965

HCPD

plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$EI.ageslope
## S = 409170, p-value = 9.57e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4534636

Spatial permutation (spin) based p-value

EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$myelin.annualroc
## S = 1071660, p-value = 7.233e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4314374

Spatial permutation (spin) based p-value

myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$amplitude.agedecline
## S = 378150, p-value = 7.236e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4136816

Spatial permutation (spin) based p-value

amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.hcpd, myelin.thalamus.hcpd, amplitude.thalamus.hcpd)
plasticity.ps
## [1] 0.00010 0.00195 0.00890
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.000300 0.002925 0.008900